Local elasticity of strained DNA studied by all-atom simulations 
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Genomic DNA is constantly subjected to various mechanical stresses arising from its biological 
functions and cell packaging. If the local mechanical properties of DNA change under torsional and 
tensional stress, the activity of DN A-modifying proteins and transcription factors can be affected and 
regulated allosterically. To check this possibility, appropriate steady forces and torques were applied 
in the course of all-atom molecular dynamics simulations of DNA with AT- and GC-alternating 
sequences. It is found that the stretching rigidity grows with tension as well as twisting. The torsional 
rigidity is not affected by stretching, but it varies with twisting very strongly, and differently for the 
two sequences. Surprisingly, for AT-alternating DNA it passes through a minimum with the average 
twist close to the experimental value in solution. For this fragment, but not for the GC-alternating 
sequence, the bending rigidity noticeably changes with both twisting and stretching. The results 
have important biological implications and shed light upon earlier experimental observations. 

PACS numbers: 87.14.gk 87.15.H- 87.15. ap 87.15. ak 
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Introduction 

Internal mechanical stress is ubiquitous in the biolog- 
ically active state of double helical DNA. In eucaryotic 
cells, DNA is densely packed in chromosomes and forced 
to bend, twist and stretch by numerous protein factors 
involved in genome regulation [lj, |2|. In procaryotes, 
DNA is subjected to a constitutive unwinding torque 
maintained by special enzymes, which leads to superced- 
ing, as in a long rope with bending and twisting elastic- 
ity @, 13] • The supercoiling and, more generally, stress- 
induced DNA forms are key factors in a variety of cellular 
processes @ . For instance, the degree of supercoiling in 
bacteria changes systematically during the cell cycle and 
in response to environmental conditions, which is accom- 
panied by activation or suppression of certain genes [6| . 
The promoter sensitivity to supercoiling stems from the 
recognition of short promoter elements by RNA poly- 
merase Q. Detailed studies indicate that it probably 
does not require DNA melting nor transitions to alterna- 
tive forms [6] . In E. coli, relaxation of the superhelical 
stress simultaneously alters the activity of 306 genes (7% 
of the genome), with 106 genes activated and others deac- 
tivated |8| . The genes concerned are functionally diverse 
and widely dispersed throughout the chromosome, and 
the effect is dose-dependent. 

The physical mechanisms of such effects are under- 
stood only partially. Long DNA is well described by the 
coarse-grained worm- like chain (WLC) model [9,[l(| sup- 
plemented with harmonic twisting and stretching elas- 
ticity [Dp - tla ]. This model nicely explains the stress- 
modulated probability of looping, wrapping around pro- 
teins, and juxtaposition of distant protein binding sites 
17]. However, it cannot account for the promoter sen- 
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sitivity to supercoiling, for instance, because in this and 
many other cases the gene regulation has a strong lo- 
cal character and is dominated by sequence effects. A 
long-discussed hypothesis is that the stress may act as 
an allosteric factor in protein-DNA recognition [lj§, [l9j . 
The supercoiling arguably changes the local properties of 
DNA, as there are small proteins with single short bind- 
ing sites that can distinguish stressed and relaxed DNA 
forms [20| ; however, it is never clear what exactly is rec- 
ognized. The supercoiling torque is distributed between 
twisting and writhing so that the untwisting of the double 
helix is estimated as 1-2% [2l[, which is below the ther- 
mal noise and too small for reliable recognition. Alterna- 
tively, the action of the torsional stress may be conveyed 
through a property other than the structure of the double 
helix. For instance, the untwisting may change the elas- 
tic parameters of DNA J22T-[2J|. The supercoiled DNA is 
governed by the interplay between the local bending and 
twisting fluctuations. If the bending flexibility or the tor- 
sional stiffness vary, parameters of thermal fluctuations 
of short DNA stretches involved in recognition could be 
noticeably affected even at low levels of stress. 

The foregoing hypothesis implies that even with small 
deformations the DNA elasticity is not exactly harmonic. 
This possibility was earlier considered in relation to spe- 
cific experiments and also to explain the discrepancies in 
twisting rigidity of DNA evaluated by different methods 
|22T - r25| . Notably, it was suggested that the stretching 
forces applied in single molecule measurements and the 
bending involved in DNA cyclization can increase the ap- 
parent twisting rigidity of DNA (25J|. The DNA double 
helix tends to overwind with small stretching |26l . |27J, 
but it is not clear if bending and/or stretching affect 
the twisting elasticity. The mechanical coupling between 
deformations of different types may be very important 
for regulation. However, the most interesting for biol- 
ogy is not the overall elasticity, but the behavior of short 
specific sequences within polymer DNA. To the present, 



all experimental studies have probed only the average 
properties of long DNA, with a few reports on sequence 
effects [23, [29j and the influence of supercoiling stress 
[22H24j . For the free relaxed double helix a good conver- 
gence of the results of different experiments is obtained 
for the bending rigidity [29|, [30( . The torsional rigidity 
has been measured by multiple different techniques, but 
the results remain controversial 25|. Also, a few esti- 



mates of the stretching stiffness have been obtained from 
nanomechanics experiments with single DNA molecules 

ei sail. 



Although the local sequence-dependent DNA elastic- 
ity and possible stress effects are difficult to reveal ex- 
perimentally, they can be probed by computer simula- 
tions. All-atom molecular dynamics (MD) simulations is 
a powerful instrument particularly suitable for this pur- 
pose. Continuous improvement of forceficlds [33- 35] and 
simulation techniques [3 a . |37| have now made possible 
free MD simulations that reproduce conformational en- 
sembles of DNA in good agreement with experimental 
data [38, 39]. Calculated statistics of fluctuations in short 
DNA qualitatively agree with the WLC theory [iollilj]. 
and the values of the elastic parameters can be measured 
with good accuracy [42|, |43[ . DNA deformation is a clas- 
sical subject of molecular mechanics [44]. In several ear- 
lier investigations, all-atom MD simulations were used 
for studying deformed DNA states obtained by exter- 
nal stretching |45l - l47l ]. twisting |48l - l5Cj |. or bending [51j . 
The required deformations were produced by either po- 
tential restraints or periodical boundary constraints. A 
promising alternative method [52j applies steady forces 
and torques to short stretches of DNA. In contrast to the 
earlier approaches, this method makes it possible to eval- 
uate elastic parameters under different types and mag- 
nitudes of external stress corresponding to physiological 
conditions. This method captures linear elastic responses 
as well as the twist-stretch coupling effect under small 
torques corresponding to a physiological degree of super- 
coiling 52] . With such approaches it has been found that, 
depending upon the base pair sequence, small twisting 
torques corresponding to physiological superhelical den- 
sity can significantly change the torsional stiffness of the 
DNA double helix [53J ]. 

In this article we present the results of the first system- 
atic study of the influence of external mechanical stress 
upon the local stretching, twisting, and bending elastic- 
ity of the double helical DNA. The numerical algorithms 
described and tested in the recent reports (52|, [53| could 
be drastically accelerated through parallelization, which 
made such computations more affordable. Two double 
helical fragments were considered, with AT- and GC- 
alternating sequences, respectively. We found that the 
apparent stretching rigidity of DNA strongly depends 
upon the method used for measuring the molecule length. 
When it is obtained b y s umming base-pair steps as in 
earlier studies [40, |43|, [5J] the sign of the twist-stretch 
coupling effect appears opposite to that measured exper- 
imentally. In contrast, much better agreement with ex- 



perimental data is obtained when the length is measured 
directly via the end-to-end distance of one helical turn. 
We argue that only the latter value corresponds to the 
experimental observable. The change in the stretching 
rigidity of DNA with external stress is qualitatively simi- 
lar for the two sequences. It grows with stretching as well 
as with increased twisting. The torsional rigidity is essen- 
tially unaffected by stretching, but it varies with twist- 
ing very strongly, and differently for the two sequences. 
Surprisingly, for the AT-alternating sequence, it passes 
through a minimum with the average twist close to the 
experimental value in solution. For this fragment, but 
not for the GC-alternating sequence, the bending rigid- 
ity noticeably changes with both twisting and stretch- 
ing. The results shed light upon the earlier experimen- 
tal observations [22h24| and have important implications 
for the possible mechanisms of allosteric gene regulation 
0111. 



Methods 

Simulation protocols 

Tetradecamer DNA fragments were modeled with AT- 
alternating and GC-alternating sequences. A dodecamer 
fragment is necessary for a full helical turn of a random- 
sequence B-DNA. The length of 14 base pairs (bp) is 
minimal for modelling of a helical turn within a longer 
DNA. This choice of the fragment length and sequences 
is consistent with and dictated by the results of the ear- 
lier studies [42], |43| . Steady stress loads were applied as 
described elsewhere [52|. This method distributes forces 
over selected groups of atoms and compensates them by 
reactions applied to other atoms so as to zero the total 
external force and torque. Because the forces are applied 
at different points internal stress and deformations are in- 
troduced that correspond to overall twisting or stretch- 
ing. The method was thoroughly verified in Brownian 
dynamics simulations of calibrated discrete WLC models 

The ranges of forces and torques are selected to com- 
prise the values used in single molecule manipulation ex- 
periments as well as the corresponding estimates for liv- 
ing cells. It is known that B-DNA becomes unstable in 
vitro with stretching forces beyond 50 pN [31;, 56]. The 
covalent bonds in long DNA are broken already with 
forces beyond 300 pN [57j |. and in living cells DNA is 
often fragmented during replication in so-called fragile 
sites (58| . A stretching load of a few tens of piconewtons 
can be exerted by a sing le molecule of RNA polymerase 
during transcription [59, 60], and forces in the nN range 
pull the chromatids during cell division [61]. The range of 
torques that do not destroy B-DNA in single molecule ex- 
periments is from -10 to +35 pNnm [621 ] . The lower limit 
is close to the estimated torsional stress due to natural 
negative supercoiling in procaryotes. These data con- 
cern the integral stability of long random sequence DNA. 



Short stretches of B-DNA can tolerate much stronger tor- 
sional strain. For instance, the DNA twisting observed in 
complexes with some bacteriophage repressors [63j cor- 
responds to torques beyond 100 pNnm. 

The classical MD simulations were carried out by run- 
ning independent trajectories in parallel on different pro- 
cessors for identical conditions. The number of processors 
varied between 32 and 48. Trajectories with the lowest 
loads started from the final states of free dynamics. The 
amplitudes of forces and torques were increased gradually 
so that simulations with higher values started from the 
final states obtained under the preceding lower values. 
The initial 0.5 ns of every sub-trajectory were discarded, 
which was sufficient for re-equilibration. 

The AMBER98 forcefield parameters^,!!!] were used 
with the rigid TIP3P water model 65]. The electro- 
static interactions were treated by the SPME method 
[37j ]. To increase the time step, MD simulations were 
carried out by the internal coordinate method (ICMD) 
[6a . l67j , with the internal DNA mobility limited to essen- 
tial degrees of freedom. The rotation of water molecules 
and DNA groups including only hydrogen atoms were 
slowed down by weighting of the corresponding inertia 
tensors [68|, |69| ■ The double-helical DNA was modeled 
with all backbone torsions, free bond angles in the sugar 
rings, and rigid bases and phosphate groups. The ef- 
fect of these constraints is insignificant, as was previ- 
ously checked through comparisons with standard Carte- 
sian dynamics [4Ct |68[. The time step was 0.01 ps and 
the DNA structures were saved every 5 ps. All trajecto- 
ries were continued to obtain the sampling corresponding 
to 164 ns of continuous dynamics, that is 2 15 points for 
every value of force (torque) . 

Additional technical details including preparation of 
initial states, treatment of rare events, evaluation of sta- 
tistical errors, and others are described elsewhere [70]. 



Evaluation of elastic parameters 

The DNA elasticity is conveniently characterized by 
three persistence lengths (PLs) corresponding to bend- 
ing, twisting and stretching that we denote here as lb, h, 
and l s , respectively. These parameters can be extracted 
from simulated canonical conformational ensembles by 
using the WLC theory that provides linear relationships 
of the following form 



D X (L) 



L 



(1) 



where L is the DNA length and x stands for b, t, or s. The 
WLC deviations D X {L) are computed from appropriate 
canonical averages as 



D b (L) = - In ((cos [0 (It)])) 
D t (L) = D [*(£)] 
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where 9(L) and &(L) are the angles of bending and twist- 
ing, respectively. The angular brackets denote the canon- 
ical averaging and D with square brackets refer to the 
variance of the variable in the brackets. The sampled con- 
formations of the double helix were analyzed by the pro- 
gram 3DNA [71(. Because the elastic parameters should 
be preferably estimated by using integral numbers of he- 
lical turns :43], only 11 central base pair steps (bps) were 
considered (central dodecamers referred to as TA 6 and 
CG6, respectively). In the following text, symbols 9, $, 
and L denote the corresponding parameters of one helical 
turn. 

According to the standard convention [72J], every base 
pair is characterized by a local Cartesian frame, with 
the xy-plane parallel to the base pair and z-vectors di- 
rected along the DNA. The bend angle 9 is measured 
between the z-vectors constructed at the opposite ends 
of a helical turn. Earlier it was shown that this mea- 
sure of bending is adequate for integral numbers of he- 
lical turns 43]. The torsional fluctuations were probed 
by three alternative methods. The end-to-end twist, $', 
was evaluated similarly to the local twist |73| . but using 
the two terminal reference frames. The cumulated local 
twist <fr" is obtained by summing the local twist at all 
base-pair steps. The last angle $'" is computed similarly 
by using the base-pair twist with respect to the optimal 
helical axis. The fluctuations of the DNA length were 
also evaluated by using three alternative methods. The 
end-to-end distance, L' was measured directly between 
the origins of the terminal reference frames. The contour 
length L" was measured by summing the distances be- 
tween the consecutive frames. The last value, L'" was 
obtained by summing the local rise from the 3DNA out- 
put. These three methods give different average L values 
and it is not evident which of them is the best estimate 
of the macroscopic DNA length. Therefore, in Eq. (p} 
we used L computed as 11 • 0.335 nm, that is by using 
the experimental length for one bps. This can cause a 
systematic bias in the measured PLs, but does not affect 
qualitative trends. 



Results 



Two stretching rigidities of the double helix 

The length of the double helix is usually evaluated 
by summing the helical rise along the molecule [4CJ, l54j . 
The rise can be measured with respect to the helical axis 
(global rise) or between the base pair frames (local rise). 
In both cases it is sensitive to alg orithmic differences be- 
tween the analysis programs 74] , and the corresponding 
l s values sometimes diverge very significantly 43]. To 
get reliable estimates we tested several possibilities and 
three representative techniques outlined in Methods are 
compared below. The end-to-end distance, L' , is a di- 
rect measure that is adequate in our case because for 
very short DNA the length fluctuations are dominated 



TABLE I Reference zero stress values of the main parame- 
ters discussed in the text. The DNA length, L, the twist, 
<3>, and the corresponding PLs, l B and It, respectively, were 
measured by three alternative methods outlined in Methods. 
The statistical errors for L and $ were about 0.4 deg and 
0.04 A, respectively. For l s and It the relative errors were 
about 4% and 5%, respectively. 



L(A) 


£ s (nm) 


$ (deg) 


k (nm) 


method TA 6 CG 6 


TA 6 CG 6 


method TA 6 CG 6 


TA 6 CG 6 


L' 34.2 35.4 
L" 38.7 38.2 
L'" 36.3 36.4 


78 172 
230 238 
342 394 


$' -5.4 15.8 
$" 349.4 372.6 
$"' 365.4 384.0 


121 123 
102 116 
140 122 



by stretching [75|, l76|. The second parameter, L", is the 
length of the three-dimensional zigzag line through the 
origins of the reference frames. By construction, L 1 < L" 
(see Table IJ). The cumulated local rise, L"', was used 
in the earlier studies [43, [54|. The local rise is one of 
the orthogonal projections of the distance between the 
neighbor frames, therefore, V" < L" . A similar value 
computed with the global rise is not considered here. 

Fig. [TJshows the extension-vs-force plots obtained with 
the above three lengths. All three plots are approxi- 
mately linear, in good agreement with the harmonic ap- 
proximation, but the V value grows much faster than 
the other two. Only the vertical positions of the theo- 
retical straight lines were fitted to the data points while 
the slopes were computed independently, which gives an 
additional check of self-consistency. The increase of L" 
is similar to that of L'" notwithstanding the divergence 
of their absolute values (see Table |T|, and this increase 
agrees with the l s value obtained from equilibrium fluc- 
tuations of L'" rather than L". To explain these observa- 
tions, note that the zigzag probed by L" forms a helical 
trace that winds around the straight segment measured 
by V '. Fig. Q] suggests that the strokes of the zigzag can 
be considered inextensible, and the end-to-end distance 
grows mainly due to flattening of angles. The local heli- 
cal parameters are obtained by decomposing each stroke 
of the zigzag into rise, shift and slide [731. All three of 
them contribute to the fluctuations of L", however, only 
the rise is affected by the applied force because the other 
two correspond to displacements nearly orthogonal to the 
force. This explains why the growth of L" is better de- 
scribed by l s obtained from L'". 

According to Fig. [TJthe double helix is characterized by 
two qualitatively different stretching rigidities. Parame- 
ter l s corresponding to fluctuations of L' can be mea- 
sured experimentally. In the experimental literature the 
stretching stiffness is conventionally characterized by the 
modulus Yf related to l s as 



I, 



Y f /3.4nmV 




The experimental estimates of Yf are around 1100 pN 

14, 31, 32], which corresponds to l s =78 nm, in reasonable 

agreement with l s computed from MD data (see Table HJ. 



10 20 30 40 50 10 20 30 40 50 
Force (pN) 

FIG. 1 Color online. DNA extension under steady stretching 
load. The molecule length was measured by three different 
methods outlined in the text. The data for L' , L" and L'" 
are shown by black dots, open red squares, and open blue 
circles, respectively. The error bars are small and they merge 
with symbols. The theoretical harmonic dependences (col- 
ored dashed lines) were plotted for the zero-stress values of 
stretching PL presented in Table U with the vertical shifts 
fitted to the data points. The left and right panels exhibit 
the results for TAg and CG6, respectively. 



This conclusion is corroborated by Fig. [2]that shows how 
the measured DNA length changes with forced twisting. 
According to experiments [26|, l27[ small twisting should 
cause extension of the double helix. It is seen that the 
end-to-end length L' indeed grows with small twisting in 
quantitative agreement with the experimental estimate. 

1.2 
.< 0.8 
^ 0.4 

-!■ o.o 

-0.4 

-40 40 80 -80 -40 40 
Torque (pNnm) 

FIG. 2 Color online. Variations of the DNA length under 
steady twisting torques. The three plots in each panel cor- 
respond to alternative definitions of DNA length explained 
in the text. The notations are same as in Fig. [T] The dotted 
traces show the expected dependences for a harmonic twist- 
stretch coupling with parameters measured experimentally 

El- 

In contrast, the stretching rigidity characterized by 
parameter l s is similar for both sequences, but signifi- 
cantly larger than the experimental estimate (see Table 
|T|. Thermal fluctuations of the local rise involve pertur- 
bations of base-pair stacking, therefore, Z s specifically 
characterizes the strength of stacking interactions. How- 
ever, this stretching rigidity is not probed in experiments. 
Fig. [2] reveals that the lengths measured by parameters 
L" and L'" both decrease with twisting, in qualitative di- 
vergence from L 1 and experimental observations [26i [27JJ . 

The stretching rigidity does not remain constant with 
forced stretching and twisting. Fig. [3] reveals that in 
stretched DNA, l a and l s deviate in opposite senses. 
The l s value corresponding to experimental measure- 
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FIG. 3 Color online. Effect of tension upon the stretching 
rigidity of DNA. The stretching PL was evaluated from the 
length fluctuations measured by three alternative methods 
outlined in the text. The reference zero stress values are 
presented in Table [I] The notations correspond to Fig. [T] 
namely, the data for l s , l s , and l s are shown by black dots, 
open red squares, and open blue circles, respectively. 



merits grows. Therefore, the molecule should gradually 
become stiffer until the stretching force approaches the 
limit of about 70 pN where the B-DNA is known to loose 
stability [73]. The stiffening agrees with the deviations 
of black points in Fig. Q] from the linear plots corre- 
sponding to the harmonic approximation. Mechanisti- 
cally, the growth of l g can be rationalized by noting that, 
with the zigzag angles flattened, the end-to-end distance 
L' approaches the zigzag length L". Since L' can never 
exceed L", the fluctuations of V should decrease, that 
is, l s grows approaching l s . The simultaneous decrease 
of l s reflects gradual weakening of base-pair stacking. 
Twisting also increases the stretching rigidity (see Fig. 
[3J. However, untwisting of TAg changes l s only slightly, 
suggesting that it passes through a minimum with torque 
r=-20 pNnm. Interestingly, the value of l s reached with 
untwisting of CGg is similar to that of TA@ . 
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FIG. 4 Color online. Effect of forced twisting upon the 
stretching rigidity of DNA. The three plots in each panel cor- 
respond to alternative definitions of DNA length explained 
in the text. The notations are same as in [3] 



Torsional rigidity 

In the previous report [53[ the torsional rigidity was 
evaluated by using the twist angle <&'" (see Methods). 
This parameter depends upon the construction of an op- 
timal straight helical axis, which can add a spurious noise 



due to bending deformations of the double helix. For 
verification, here the torsional dynamics are analyzed by 
three alternative methods including the earlier one. The 
end-to-end twist $' is most appropriate for comparisons 
with experiment because it closely corresponds to that 
measured in experiments with long DNA. The cumulated 
local twist $" represents another reasonable alternative 
and it was added as an additional check. 
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FIG. 5 Color online. DNA twisting under steady torques. 
The twist angle was measured by three different methods 
outlined in the text. The data for $', <&", and <&'" arc 
shown by black dots, open red squares, and open blue cir- 
cles, respectively. The error bars are small and they merge 
with symbols. The theoretical harmonic dependences (col- 
ored dashed lines) were plotted for the corresponding zero- 
stress values of the torsional PL presented in Table U with 
the vertical shifts fitted to the data points. The left and right 
panels exhibit the results for TAg and CGg, respectively. 
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FIG. 6 Color online. Effect of stretching upon DNA twist- 
ing. The three plots in each panel correspond to alternative 
measurements of DNA twisting as explained in the text. The 
notations are same as in Fig. \E\ The dotted traces show 
the expected dependences for a harmonic twist-stretch cou- 
pling with parameters measured experimentally for random- 
sequence DNA 26] . 



The external torque changes DNA twisting as shown 
in Fig. [5] In contrast to stretching, the three alterna- 
tive measures of angle $ give very similar results in spite 
of the divergence of the reference zero-stress values (see 
Table HJ. Similarly to Fig. [lj only the vertical positions 
of the theoretical straight lines were fitted to the data 
points, with the slopes computed independently This ad- 
ditionally checks the self-consistency and one may note 
that the deviations from the harmonic law are smaller 
for $' and $'" than for $". Earlier single-molecule ex- 
periments revealed that DNA overwinds when stretched 



[2 a . |27| | . This effect is well reproduced with any of the 
three methods (see Fig. [5]) sometimes with good quan- 
titative agreement. The dashed lines in Fig. [6] represent 
the experimental dependence for small forces below 30 
pN [26j. With stronger extension the twist should start 
to fall. For TA6 this experimental observation is not re- 
produced, but for CGg a transition from an ascending 
trend to an irregular plateau is indeed observed at about 
30 pN. This irregular dependence is not due to errors or 
hidden statistical noise. For verification, we reduced the 
force from 50 to 40 pN, repeated the MD simulations, 
then raised the force back to 50 pN, and carried out one 
more run. The results of this back and forth test were 
within the error limits shown in Fig. [6] 
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FIG. 7 Color online. Effect of forced twisting upon the tor- 
sional rigidity of DNA. The torsional PL was evaluated from 
the twist fluctuations measured by three alternative meth- 
ods outlined in the text. The reference zero stress values are 
presented in Table [TJ The notations correspond to Fig. [5] 
namely, the data for l t , l t , and l t are shown by black dots, 
open red squares, and open blue circles, respectively. 



The measured torsional rigidity changes with forced 
twisting as shown in Fig. [7] The three alternative mea- 
sures of twist yield very similar results all showing strong 
variations of l t , with a remarkable qualitative difference 
between the two sequences. These rigidity profiles agree 
with the non-linear features of the $(r) plots in Fig. \5\ 
Indeed, for CG6 they are concave and for TA 6 the har- 
monic law corresponding to the zero-stress rigidity over- 
estimates the twisting of both signs. The twisting rigid- 
ity of CGg grows steadily in the whole range of torques 
tested. In contrast, for TA 6 an opposite trend is observed 
under small torques, but l t passes via a minimum under 
positive torques. A qualitatively similar behavior was 
experimentally observed for one natural DNA sequence 

mm. 

The growth of rigidity with torques of both signs agrees 
with simple physical intuition for a twist energy profile 
resembling a flat-bottomed basin with vertical walls. In 
this case the system cannot go very far even with strong 
energy fluctuations. The range of torques applied to 
CG6 was extended to check the existence of a minimum 
under negative torques. It is seen, however, that the 
minimum is not reached although the decrease of It be- 
comes less steep with untwisting. This behavior indicates 
that anomalously frequent strong untwisting fluctuations 



should occur in GC-alternating DNA under normal tem- 
perature. 
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FIG. 8 Color online. Twisting rigidity of DNA under small 
and large torques. The results of the earlier report 53] (open 
symbols) are compared with those of the present study (filled 
symbols). Squares and circles show the data data for TAg 
and CG6, respectively. The DNA twisting was evaluated by 
using angle <&"'. 

The results in Fig. [7] confirm and corroborate the con- 
clusions of the previous report where smaller torques were 
considered [53j . The earlier data are compared with those 
of the present study in Fig. [H It is seen that the two se- 
ries of simulations are consistent in spite of the differences 
in protocols. Each open circle and open square in Fig. 
H] correspond to a single continuous trajectory, therefore, 
this figure confirms ergodicity and validates the much 
faster protocol introduced here. The new plots also look 
less noisy, which can be attributed to the absence of slow 
non-canonical a/7 dynamics. In the previous calcula- 
tions, such transitions occurred almost exclusively in ter- 
minal bps [52, [53|, nevertheless, they affected the middle 
fragments allosterically and contaminated the results. 

The strong torsional anharmonicity is not seen in the 
shapes of the probability distributions of twisting fluctu- 
ations of the whole fragment. These distributions remain 
nearly Gaussian, with the widths changing in agreement 
with Fig. [JJ(see Refs. [53,(70). In contrast, the pattern of 
single-step twist fluctuations qualitatively explains the ef- 
fect revealed in Fig. [JJandUl As seen in Fig. [H with a no- 
table exception of the adenine-phosphate-thymine steps 
(ApT), these distributions strongly differ from Gaussians 
predicted for harmonic WLC model with the measured 
It values. Surprisingly, for TpA and CpG steps these 
shapes qualitatively change with twisting. With nega- 
tive torques the distributions in the upper two panels are 
strongly positively skewed, but they gradually become 
negatively skewed as the torque changes the sign. The 
same is true for the GpC distributions although in this 
case the effect is much smaller. Some of the TpA and 
CpG distributions exhibit clear humps suggesting that 
the twisting in these steps is best described by double- 
well potentials with low transition barriers. 

The CpG and GpC distributions in Fig. [9]behave simi- 
larly, that is, they become wider with untwisting in qual- 
itative agreement with the lt(r) plots for CGg in Fig. [5] 
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FIG. 9 Color online. The normalized probability densities 
of twisting fluctuations in the four types of bps represented 
in TA6 and CG6 fragments. Open circles show the com- 
puted distributions for strong negative (left, red), medium 
(middle, green), and strong positive (right, blue) values of 
external torque, respectively. The corresponding torque val- 
ues were r= -40, +40, and +80 pNnm for TAg and t= -40, 
0, and +40 pNnm for CG6- The black solid lines show the 
analytical Gaussian distributions for It measured under the 
middle torque values. These curves are placed to match the 
maxima of the corresponding computed distributions. 



In contrast, in TAg the two alternating dinuclcotidc steps 
behave differently. The width of the ApT step distribu- 
tions changes monotonously in the whole range of torques 
probed, that is, the minimum of It at 40 pNnm in Fig. [5] 
is exclusively due to TpA steps. The ApT distributions 
also exhibit a striking feature. The centers of all plots are 
shifted in agreement with the sign of the applied torque, 
however, the magnitude of the shift is small compared 
to the change in the distribution width. As a result, the 
probabilities of strong untwisting fluctuations are higher 
with t=+80 pNnm than with r=-40 pNnm. The effect 
is small, but statistically significant (see also Ref. [70h . 
This feature is counterintuitive because it cannot be re- 
produced with the WLC model. 
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FIG. 10 Effect of stretching forces upon the torsional rigidity 

of DNA. The measured values of l t are shown for TAg and 

CG6 on the left and right panels, respectively. 

In contrast to twisting, small stretching has virtually 



no effect upon the torsional rigidity of DNA. The corre- 
sponding data are shown in Fig. 1101 For clarity, only 
the l t values are shown. The variations are small and 
rarely exceed the statistical errors. When l t is measured 
by using magnetic tweezers the common stretching load 
is smaller than 20 pN [62[ and the data in Fig. Qjj] indi- 
cate that it can noticeably affect the results only due to 
mechanisms that are not reproduced in the present DNA 
model. 



Bending rigidity 

In long DNA, stretching naturally flattens bends, 
whereas twisting causes looping and supercoiling, that 
is, increases bending in some DNA stretches. These ef- 
fects are strong; the accompanying changes in the bend- 
ing rigidity are hardly measurable experimentally and 
this possibility usually is not considered. The atom-level 
modeling is the only currently available method that can 
check whether or not the bending rigidity of DNA in 
principle can be affected by the twisting and/or tensional 
stress. The results of the first such tests are shown in Fig. 
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FIG. 11 Color online. Variation of the apparent bending PL 
with forced stretching and twisting. The results for CG6 and 
TAg are shown by red squares and blue circles, respectively. 



The measured bending PL of CGg exhibits only small 
variations with both twisting and stretching. In contrast, 
for TAq these variations significantly exceed statistical 
uncertainty and reveal interesting trends. Notably, the 
right panel of Fig. [TTJ reveals that bending in TAg in- 
creases with stretching, which is opposite to the expected 
flattening effect. At the same time, the lb(r) dependence 
in the left panel of Fig. [TTJ passes through a maximum at 
40 pNnm, that is exactly where the torsional PL reaches 
the local minimum in Fig. [7J A closer look reveals 
that these trends are accompanied by subtle qualitative 
changes in the bending dynamics. By using the base pair 
coordinate frames provided by the program 3DNA [71( 
one can conveniently characterize the bend direction as 
follows. Consider two coordinate frames constructed at 
the first and the last base pairs, respectively. According 
to the standard convention [72[ , the two xy-planes dissect 
the double helix approximately parallel to the base pair 
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planes. The corresponding two z-vectors approximate the 
local directions of the helical axis. If the z-vectors are 
not parallel we can construct the orthogonal projection 
of the second z- vector upon the first xy-plane. The spher- 
ical azimuth angle ip is measured between the projected 
z-vector and the x-vector of the projection plane. With 
the x-vector corresponding to the standard convention 
J72] , the value of <p is close to zero when the molecule is 
bent towards the minor groove in the middle of the he- 
lical turn. A few representative distributions of angle tp 
are shown in Fig. Q21 
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FIG. 12 Color online. Effects of stretching and twisting upon 
the anisotropy of bending. The normalized probability dis- 
tributions for the bending azimuth angle ip are shown for 
TAg twisting (left panel) and stretching (middle panel), and 
for CG6 twisting (right panel). The red, blue, and green 
dots correspond to torques r= -40, +40, and +80 pNnm 
(left panel), forces of 0, 20, and 50 pN (middle panel), and 
torques r= -80, 0, and +60 pNnm (right panel). 

The red distribution in the middle panel indicates that 
the unstressed TAg fragment prefers to bend towards the 
minor groove. The origin of this anisotropy should be 
studied additionally because it is probably inherent in the 
overall dynamics rather than caused by local end effects 
or construction of the coordinate frames. Here we use it 
just as an indicator. It is seen that the original anisotropy 
increases with both stretching and unwinding, but pos- 
itive torques reduce it. As a result, with r=40 pNnm 
the azimuth distribution becomes even, and with further 
twisting the anisotropy of an opposite sign appears (left 
panel). This behavior is in remarkable contrast to that 
of CG6- For CG6 the bending is also preferable towards 
the minor groove, but twisting only causes rotation of 
this direction in agreement with the relative orientation 
of the minor groove in the middle of the fragment (right 
panel). 

Comparison of Fig. [T^] with Fig. [7J and O suggests 
that r=40 pNnm corresponds to a transition state be- 
tween two qualitatively different dynamic patterns and 
that this behavior is attributable to the specific proper- 
ties of TpA steps. The results in Fig. [7] and M can be 
readily rationalized and qualitatively reproduced in an 
appropriate coarse-grained model with local twisting de- 
scribed by a double- well potential. Global bending of the 
double helix results from local deviations of bps geome- 
try described by parameters Roll, Tilt, Slide, and Shift 



72]. Analysis shows that, in TpA steps, all of them are 
affected by twisting. Roll and Slide change more than 
other, with Slide exhibiting a bi-modal pattern of fluctu- 
ations [70(. These data demonstrate that a mechanical 
link between twisting and bending is an inherent prop- 
erty of TpA dinucleotides, therefore, we can qualitatively 
explain the results in the left panel of Fig. Q21 The 
quantitative relation is much more difficult to establish 
because global bending results from a complex summa- 
tion of local motions over the whole fragment, including 
correlations and helical rotation. The correspondence of 
the transition states for twisting and bending in TAg may 
be a coincidence, nevertheless, the results in Fig. [TTIand 
H2] evidence that twisting and stretching can produce un- 
expected sequence dependent effects upon the bending 
dynamics in DNA. The stress response is complex and it 
cannot be reduced to an altered bending rigidity of the 
underlying WLC model. 



Discussion 

A very good agreement of experiments on polymer 
DNA with the WLC model flMl(fy have led to an ex- 
aggerated belief that the harmonic approximation is suf- 
ficient for describing all essential properties of the DNA 
double helix. In fact, these remarkable results cannot be 
considered as evidences of harmonicity because the ad- 
ditive ladder construction of the double helix effectively 
hides local heterogeneity and anharmonicity. Due to this 
additivity, and the central limiting theorem of the prob- 
ability theory, various experimental data converge to the 
WLC model regardless of the local DNA properties, with 
only a few concatenated bps being sufficient for the ap- 
parent statistical equivalence with a harmonic elastic rod 
[53|. This effect shadows the true mechanical properties 
of the DNA double helix which remain elusive. 

The present study evidences that, under normal tem- 
perature, the local DNA elasticity is strongly anhar- 
monic, in agreement with the early hypotheses [78| and 
some experimental data [23 - t24J j. The results of compu- 
tations using empirical forcefields certainly require fur- 
ther verification. New experimental approaches need to 
be developed for this purpose because currently available 
methods can probe only the average elastic parameters 
of long molecules. 

The computed values of all elastic parameters reason- 
ably agree with the data for polymer DNA obtained by 
different experimental methods. The earlier controversy 
concerning the stretching (Young's) modulus [43[ is clari- 
fied here by comparing different procedures for measuring 
the length of the double helix. The experimental bending 
rigidity is characterized by lb ~ 50 nm [30j. The mea- 
sured It values vary between 36 and 109 nm depending 
upon specific methods and conditions [25j , The stretch- 
ing PL is about 80 nm 0, [M H3- MD simulations 
give somewhat larger values, that is, the DNA stiffness 
is slightly overestimated. [40, Hj] This discrepancy is 



not large and it can be attributed to a combination of 
factors like inexact correspondence between the micro- 
scopic geometric parameters and experimental observ- 
ables, the neutralizing salt conditions in MD, and the 
small size of the modeled fragments predictably leading 
to strong sequence and end effects. As shown here, MD 
also quantitatively reproduce the reciprocal coupling be- 
tween twisting and stretc hing revealed in recent magnetic 
tweezer experiments [26l . |27| . The overall agreement is 
quite remarkable because none of the MD forcefield pa- 
rameters was adjusted to fit the computed DNA elasticity 
to experiment. One may reasonably hope, therefore, that 
the detailed microscopic picture provided by simulations 
captures the qualitative physical trends dictated by the 
atom-level mechanics of the double helix. 

Our results indicate that the most significant anhar- 
monicity is inherent in the torsional DNA deformations, 
which is attributable to the special character of stack- 
ing interactions. The twisting occurs due to sliding 
within the stacks; this motion is essentially barrierless 
and its amplitude significantly exceeds the zone where 
the harmonic approximation is valid. Even small twisting 
torques can cause significant changes in elastic parame- 
ters. The qualitative difference in the stress response of 
the torsional rigidity of AT 6 and CG6 indicates that this 
property is strongly sequence-dependent. Opposite local 
trends can mutually cancel out, which makes difficult de- 
tection of anharmonic effects in long DNA. There are a 
few reports in the literature where relevant experimen- 
tal data qualitatively differ from predictions of harmonic 
models. This occurred with some natural plasmid DNA 
]24j and also with synthetic alternating sequences (18J . 
The latter were recently found anomalous as regards the 
sequence-dependent bending rigidity [29J]. These earlier 
results require additional investigations. 

The mechanical strain is an ubiquitous attribute of 
living DNA and a key factor in genome packaging and 
regulation. The common magnitudes of natural forces 
and torques are quite large, therefore, a wide spectrum 
of non-linear structural responses should be anticipated, 
with elastic deformations at one end of the scale, and lo- 
cal melting at the other end. A few anharmonic effects 
revealed here have some interesting implications for gene 
regulation mechanisms. According to Fig. HI with the 
helical twist slightly shifted from the equilibrium value 
the sequence dependence of the DNA elasticity can be 
significantly changed and enhanced. The measured tor- 
sional stiffnesses are similar without applied torque, but 
diverge with both twisting and untwisting. For other se- 
quences, similar behavior can be anticipated for bending 
and stretching. The deformability of DNA is long con- 
sidered as a possible governing factor in the sequence- 
specific site recognition [79j , but this mechanism requires 
strong sequence dependence of local elastic parameters 
compared to that observed in experiments with long free 
DNA [28]. As we see the properties of the relaxed DNA 
cannot be simply transferred to supercoiled and/or pro- 
tein bound DNA states. Additional studies are necessary 



to check if the elastic properties of the specific binding 
sites are sensitive to external stress. 

Unexpectedly, we found that the torsional rigidity of 
AT6 passes via a minimum under moderate positive 
twisting torques. This feature is probably due to a bi- 
modal character of twist fluctuations in the TpA steps 
(Fig. |H]). The average twist of AT 6 with r=40 pNnm ac- 
tually is very close to the experimental value in solution 
[80l |8l| because in free AMBER simulations the DNA 
structures are somewhat underwound 64]. In this state 
the TpA steps exhibit a distribution of twist fluctuations 
corresponding to a saddle point between two domains of 
attraction (Fig. |H|). This point also coincides with the 
maximum in the measured bending PL accompanied by 
inversion of the local bending anisotropy. 

Earlier it was suggested that the TpA steps can adopt 
at least two distinct conformational states. Depending 
on the sequence context, there is always a temperature 
range where the TpA steps exhibit slow conformational 
transitions with relaxation times beyond the nanosecond 
time range [821 ] . These slow motions should involve ex- 
tended DNA stretches, that is, these are global transi- 
tions accompanied by switching in the TpA steps. The 
same local switching is probably responsible for the un- 
usual effects observed here. The exceptional properties of 
the TpA steps are long-known in molecular biology [83l j . 
These steps are found in both narrowings and widen- 
ings of the minor B-DNA groove [8J, [85| ■ Periodically 
spaced TpA steps is the most statistically significant fea- 
ture of DNA sequences that provide optimal DNA wrap- 
ping around nucleosome particles |8g]. Switching of lo- 
cal bending anisotropy in response to variable torsional 
stress may play some role in the control of DNA wrapping 
and unwrapping. Future studies will show whether or not 
these processes are related with the unusual microscopic 
dynamics revealed in our computations. 

According to Fig. [9] the strong variation of the twist- 
ing rigidity of CGg is mainly due to CpG steps. They 
exhibit anomalously high probability of negative twist 
fluctuations with torques around zero. The CpG steps 
are found in a number of known protein binding sites, 
but their most important biological role is related with 
C5-cytosine methylation and epigenetic regulation mech- 
anisms [83] ■ The recognition of CpG sites is a complex 
multi-facet process because they exist in three methy- 
lation states with distinct functions and because spe- 
cific binding, methylation and demethylation can occur 
on both free and nucleosome bound DNA [88l - [9l| . In- 
terestingly, methylation of free DNA strongly depends 
upon supercoiling, with the superhelical density acting 
smoothly in a dose-dependent manner [92j. The corre- 
sponding catalytic mechanism requires cytosine flipping 
from the DNA stack into a protein pocket [93[ . The low 
energy pathway of this flipping transition may require a 
strong twisting fluctuation of the CpG step, which would 
explain the effect of the torsional strain [92| . 

The above specific examples suggests a more general 
hypothesis concerning the possible role of strong DNA 
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fluctuations in gene regulation, with the non-linear elas- 
ticity as the governing factor. There are many long- 
known and well-documented processes in vivo where 
strongly deformed conformations are involved instead 
of canonical B-DNA. Deformed DNA conformations are 
ubiquitous in X-ray structures of protein-DNA com- 
plexes, so that one may wonder why there is no evo- 
lutionary pressure towards proteins that can recognize 
relaxed B-DNA? It was shown that the activity of pro- 
moters regulated via strongly deformed DNA states can 
be increased by mutations that reduce the deformation 
energy [53, [SJj , but these mutations are not selected in 
vivo. It is possible that the prevalence of large DNA de- 
formations is not a trivial consequence of its flexibility, 
but a necessity of regulatory mechanisms that involve me- 
chanical stress. The larger the deformation - the lower its 
probability and the population of such state. However, 
these low probabilities can strongly change in response 
to small regulatory impulses, in contrast to populations 
of low energy states. The non-linear elastic effects should 
play an important role in such regulation because they 
can greatly amplify the input signal and also make possi- 
ble complex responses like coupling of the amplitude and 
the anisotropy of local bending to the torsional stress as 
in the TA 6 fragment studies here. Similar ideas were dis- 
cussed in the earlier literature. This hypothesis is com- 
plementary to the view of DNA as an allosteric protein 
cofactor [95] used to explain the smooth modulation of 
gene activity during cell development 96] . The effects of 
mechanical strain upon the probabilities of strong fluctu- 
ations in DNA represent significant interest and require 
further studies. New insights in this direction can be ob- 
tained by using MD simulations of DNA in steady stress 
conditions [52] and this work is continued. 
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APPENDIX 



Simulation protocols 



Tetradecamer DNA fragments were modeled with AT- 
alternating (ATATATATATATAT) and GC-alternating 
(GCGCGCGCGCGCGC) sequences. The choice of the 
fragment length and sequences is consistent with the re- 
cent computations [42, Sj, [52| and it was dictated by the 
following considerations. An integral number of helical 
turns is preferable for evaluation of the elastic parame- 
ters of DNA and one helical turn is optimal because of the 
rapid growth of the principal relaxation times with the 
chain length [431 ]. These molecules are homopolymers of 
ApT and GpC dinucleotides, therefore, they cannot have 



distinguished asymmetric structures like static bends. 
True homopolymer DNA duplexes have special proper- 
ties and, in free MD with the AMBER forcefield, these 
structures deviate from the canonical B-DNA stronger 
than AT- and GC-alternating sequences [97j . 

The classical MD simulations were carried out by run- 
ning independent trajectories in parallel on different pro- 
cessors for identical conditions. The number of processors 
varied between 32 and 48. The starting states were pre- 
pared as follows. The solute in the canonical B-DNA con- 
formation [98] was immersed in a 6.2-nm cubic cell with 
a high water density of 1.04. The box was neutralized by 
placing Na + ions at random water positions at least 5 A 
from the solute. The system was energy minimized and 
dynamics were initiated with the Maxwell distribution of 
generalized momenta at low temperature. The system 
was next slowly heated to 293 K and equilibrated during 
1.0 ns. After that the water density was adjusted to 0.997 
by removing the necessary number of water molecules se- 
lected randomly at least 5 A from DNA and ions, and 
the simulations were continued with NVT ensemble con- 
ditions. Independent starting states for parallel trajec- 
tories were prepared by redistributing the counterions 
around DNA and re-equilibrating the system with dif- 
ferent random sets of initial momenta. The statistical 
independence was verified as explained further below. 
Steady stress loads were applied as described elsewhere 
[52j . Trajectories with the lowest loads started from the 
final states of free dynamics. The amplitudes of forces 
and torques were increased gradually so that simulations 
with higher values started from the final states obtained 
under the preceding lower values. The initial 0.5 ns of 
every sub-trajectory were discarded, which was sufficient 
for re-equilibration. 

The AMBER98 forcefield parameter sj33,l6J| were used 
with the rigid TIP3P water model 65]. The electro- 
static interactions were treated by the SPME method 
37| . with the common values of Ewald parameters, that 
is 9 A truncation for the real space sum and j3 « 0.35. 
The temperature was maintained by the Berendsen algo- 
rithm 99] applied separately to solute and solvent with a 
relaxation time of 10 ps. To increase the time step, MD 
simulations were carried out by the internal coordinate 
method (ICMD), J6(||67| with the internal DNA mobil- 
ity limited to essential degrees of freedom and rotation of 
water molecules and internal DNA groups including only 
hydrogen atoms slowed down by weighting of the cor- 
responding inertia tensors. [63, [69| The double-helical 
DNA was modeled with all backbone torsions, free bond 
angles in the sugar rings, and rigid bases and phosphate 
groups. The effect of these constraints is insignificant, as 
was previously checked through comparisons with stan- 
dard Cartesian dynamics [40J, |68[. The time step was 
0.01 ps and the DNA structures were saved every 5 ps. 
All trajectories were continued to obtain the sampling 
corresponding to 164 ns of continuous dynamics, that is 
2 15 points for every value of force (torque). Statistical 
convergence and errors were evaluated by the method of 
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block averages (see further below). 

The terminal base pairs open rather frequently during 
nanosecond time scale MD, which significantly perturbs 
the flanking DNA structure. Because this dynamics can- 
not be averaged during the accessible duration of MD tra- 
jectories, we blocked it by applying non-perturbing upper 
distance restraints as explained elsewhere [52J]. A simi- 
lar approach was used to confine the backbone dynam- 
ics to the canonical B-DNA zones. Rare non-canonical 
7 9+ — > 7' switches accompanied by a 9 ~ /a 9+ and /3*//3 ff+ 
dyn amic s complicate statistical analysis of MD trajecto- 
ries [100l | . The population of non-canonical conformers is 
overestimated with all versions of the AMBER94 force- 
field including parmbscO J35J | most used in recent years 
101]. With the parm98 modifications [64} designed to in- 
crease the average DNA twist, this problem seems to be 
less severe [52|]. However, even with a hypothetical ideal 
forcefield, such transitions should have been blocked be- 
cause they are too rare for accurate averaging. It appears 
that the following flat-bottom restraint upon 7 torsions 
solves the problem 



U~, 



^ 7 (7-2tt/3) 2 , 
0, 



7 > 2tt/3 

-7T < 7 < 27T/3 



with Ay «30 kcal/mol. This noninvasive approach does 
not perturb the dynamics in the canonical zones, which 
makes possible comparison with the previous long-time 
simulations of similar DNA fragments. 



4 

02 

Q_ 1 




4 
_3 

?1 



-180 -90 90 180 -180 -90 90 180 
a(°) p(°) 




-180 -90 90 180 -180 -90 90 180 

in s(°) 



4 4 



^». 




-180 -90 90 180 -180 -90 90 180 

E n c, n 

FIG. 1 FIG.ITJi: Distributions of backbone torsion angles in 

free MD simulations of d(ApT)7 (red circles) and d(GpC)7 

(blue circles). 



In all simulations the B-DNA conformations were well 
conserved without signs of accumulated deformations. 



Representative probability distributions of backbone tor- 
sions are shown in Fig. [1] They correspond to standard 
B-DNA dynamics (lOOJ . 



Evaluation of statistical errors 

Evaluation of errors in MD simulations is based upon 
the following assertions from the probability theory. Con- 
sider a random variable x with expectation Mi = £ and 
variance Di 
compute 



a 1 . We can take n samples of x and 



x = — (xx + x 2 + 
n 



%n) 



1 n 



x k 



fc=l 



and 



s 2 = ^£(* fc -*) 2 
fc=i 



called the sample average and variance, respectively. 
Both x and S 2 are random variables, with Mi = £ and 
MS* 2 = a 2 , i.e. x and S 2 provide unbiased estimates of 
£ and a 2 , respectively. It is also known that 



Di 



n 



and, if x is a Gaussian random variable, 

D5 2 = 2<T 4 /(n-l). 



(3) 



(4) 



Eq. ([3]) and (Q| are used for evaluation of statistical er- 
rors. 

Consider evaluation of torsional fluctuations, for in- 
stance. In this case, the random variable is the twist 
angle of one helical turn, $, with expectation M$ and 
variance D$. The torsional persistence length is com- 
puted as It = L/D$. The canonical moments M$ and 
D$ are thermodynamic limits and they are estimated 
by using, respectively, the sample average and variance 
computed over all n points saved during an MD tra- 
jectory. However, Eq. ^ and (0| cannot be applied 
straightforwardly because they are valid only for statis- 
tically independent samples, i.e. the time intervals be- 
tween the MD states must be suitably large compared to 
the torsional relaxation time. The data saving interval 
is commonly much smaller, therefore, the errorsare eval- 
uated by using the method of block averages [l02i . 1 103J ) . 
The trajectory is successively divided in n' = 2,4, ..., 2 15 
stretches (blocks) and the sample variances S' 2 are com- 
puted by using n! block averages instead of individual 
samples. When the blocks are longer than the torsional 
relaxation time the block averages are independent and 
S' 2 /n' « const — a 2 /n, where h is the effective number 
of independent samples provided by the trajectory. This 
value should be used in place of n in Eq. ([3]) and ((4]). In 
practice, it is convenient to draw the plots of 



fi(n') 



nS' [ 
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with respect to log 2 n'. When statistical independence is 
reached, such plots exhibit a plateau with Q rs n/n = t c , 
which gives the required estimate of h. Parameter r c 
is the effective correlation time measured in trajectory 
saving steps. 

Fig. [2] shows how this works for Brownian dynamics 
(BD). We consider eight trajectories from earlier pub- 
lished simulations of a discrete wormlike chain (WLC) 
model of 14-mer DNA [H, [52| . Each data set involved 
2 15 consecutive configurations saved with a 5 ps interval 
(about 164 ns in total), that is exactly as in MD simula- 
tions. With n' decreasing, the plots display emergence of 
a plateau and the growth of statistical noise. The black 
dashed line displays the results for a single eight times 
longer trajectory. It is seen that the plateau becomes 
less noisy, but its level does not change. This confirms 
the validity of the above derivations, i.e. r c is constant 
and, consequently, Eq. ([3]) holds. Fig. [5] shows that a 
reasonably accurate estimate of the plateau value can be 
obtained with n' — 2 6 or larger. On the other hand, 
the plateau emerges with the block length 10r c or larger, 
therefore, the lower estimate of the necessary total dura- 
tion of trajectory is 640r c . 



60 













l 














^£" 


.... 










"V "'-. 


'--... 



10 
log 2 (n') 



18 



FIG. 2 FIG.[2]3: Analysis of statistical errors in BD simula- 
tions by the method of block averages 102, 103]. The solid 
color lines correspond to eight similar 164 ns trajectories. 
The black dashed line displays the results for a single eight 
times longer trajectory. 



The foregoing analysis is equally valid when several 
trajectories are concatenated and considered as a sin- 
gle longer trajectory. It is only necessary that the sub- 
trajectories are statistically independent and suitably 
longer than r c . A representative example is shown in the 
upper panel of Fig. [3J We consider the torsional fluc- 
tuations of CGg fragments. The dashed black plot was 
obtained by processing a continuous 164 ns trajectory 
form our earlier report 53[ . The solid red line shows the 
analogous results for the protocol used here (32 concate- 
nated sub-trajectories). In both cases no external stress 
was applied. The plots show good convergence except for 
small n' where the statistical noise is high. The statis- 
tical dependence between the subtrajectories is analyzed 
by using time cross-correlation functions computed as 



C 13 {t) = 



Cv(0) 
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FIG. 3 FIG. [Hp: Analysis of statistical er rors in M D simu- 
lations by the method of block averages Il02l . I103H . In the 
upper panel, the solid plots were obtained by concatenat- 
ing 32 parallel MD trajectories. The red trace shows the 
results of free dynamics without external stress. Analysis of 
a continuous 164 ns trajectory with the same total volume 
of sampling is shown by the black dashed line. The lower 
panel shows the time autocorrelation function of the $ angle 
from the continuous trajectory (black dashed trace) and the 
corresponding average cross-correlation function computed 
for 32 parallel trajectories (red circles). The blue and ma- 
genta traces in the upper panel show the results for applied 
torques of +40 and -40 pN-nm, respectively. 



C ij (t) = ($ 4 (r + t)&(T)) T - ($ l )($ J ) 

(5) 



where the superscripts i and j refer to individual sub- 
trajectories. The lower panel compares the cross- 
correlation function averaged over all i ^ j with the time 
autocorrelation function from the long trajectory. The 
latter is obtained by setting i = j in Eq. (J3J) - This figure 
confirms that the subtrajectories are statistically inde- 
pendent from the very beginning. The r c f=a 40 estimated 
from the upper panel gives n=820 and the relative error 
of 4.5% in the measured It values, which is sufficient for 
our purposes. This relaxation time is approximately two 
times that estimated from the decay of autocorrelations 
by using the condition C(t) = 1/e. A similar consistency 
is observed for BD simulations (see Fig. [5] and Ref. [43h . 

The upper panel of Fig. [3] displays two additional 
traces obtained in simulations with applied external 
torques of ±40 pN-nm. The torsional elasticity of the 
CGg fragment changes with the applied torque so that 
the double helix becomes softer with untwisting. This 
gives a noticeable increase in the torsional relaxation 
time, and, accordingly, the three traces corresponding to 
positive, zero, and negative torques diverge. As a result, 
for this DNA fragment, the relative error in the measured 
It values increases with untwisting. 
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FIG. 4 FIG.HJ3: The normalized probability density P$ un- 
der different twisting. Data for TAg and CGg are exhib- 
ited in the left and right columns, respectively. For each 
fragment, the results are displayed for three representative 
torque values. For TAg, the plots correspond to r= -40, +40, 
and +80 pN-nm, from left to right, respectively. For CG6, 
the torque values were t= -40, 0, and +40 pN-nm. The an- 
alytical distributions Eq. l(6} for to the measured values of 
It and <& T are shown by solid lines. The lower panels display 
the same data in semi-logarithmic coordinates. 



Results 

In a harmonic approximation, the canonical distribu- 
tion of twisting fluctuations reads 




-2-10 1 
Slide (A) 



2 2 

Shift (A) 



FIG. 6 FIG.[6]3: The normalized probability densities of fluc- 
tuations of base-step parameters roll, tilt, shift, and slide in 
TpA dinucleotides of TAg . The definition of parameters cor- 
responds to the standard conventio n. |72| Their values were 
computed by the program 3DNA. |7l| The red, green, and 
blue plots show the data for torque values t= -40, +40, and 
+80 pN-nm. 



P$ ~ exp 



2L 



(* - *r)' 



(6) 



where $ T is the shifted equilibrium twist angle under 
torque r. The persistence length l t defines the width 
of the distribution and it must be constant if the har- 
monic approximation is valid. The computed patterns of 
fluctuations appeared qualitatively similar for all three 
measured twist angles. In Fig. @] a few representative 
distributions are exhibited for fluctuations of angle $'. 
All of the distributions are close to the analytical Gaus- 
sians defined by Eq. ©, however, their widths vary in 
agreement with the observed changes in l t . Systematic 
deviations from the Gaussian shape are noticeable only 
with the largest positive torques. 

























N 


v„ 






n n 










%■ 




° d 


w : 










"_■ -z 














• J 



15 



20 



25 



30 



35 



40 



FIG. 5 FIG. 03: The normalized probability densities of 
twisting fluctuations in the ApT steps of the TAg fragment 
shown in semi-logarithmic coordinates. The red, green, and 
blue plots show the data for external torques of -40, +40, 
and +80 pN-nm, respectively. 
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